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1 . 0  OVERVIEW 

1.1  Introduction 

The  goal  of  the  current  research  effort  is  to  develop  a  boundary  element 
method  (BEM)  analysis  to  study  the  effect  of  plasticity  on  the  growth  rate  of 
fatigue  cracks.  The  investigation  mainly  focuses  on  the  crack  tip  behavior  of 
two-  and  three-dimensional  crack  propagation  and  closure  effects  under  cyclic 
loading  for  cracks  where  the  plastic  zone  size  is  large  relative  to  the  crack 
size.  This  problem,  known  as  small  crack  or  small  flaw  problem,  stems  from 
the  rapid  growth  of  small  surface  cracks,  hence,  the  deviation  from  linear 
elastic  fracture  mechanics  (LEFM)  predictions. 

The  United  States  Air  Force  has  recognized  the  importance  of  the  small 
flaw  modeling  problem  and  Is  currently  funding  an  experimental  study  of 
turbine  engine  materials  at  SwRI.  The  analytical  study  reported  here 
supplements  this  effort  to  provide  a  most  complete  utility  of  the  work  for  the 
Air  Force. 

Because  of  many  unresolved  issues  related  to  the  two-dimensional 
analysis,  a  major  portion  of  the  current  study  is  focused  on  two-dimensional 
modeling  issues.  A  significant  issue  in  the  study  of  fatigue  crack  growth  is 
the  realization  of  crack  closure  phenomenon,  first  proposed  by  Elber  [1,2], 
where  the  residual  plastic  stretch  left  in  the  wake  of  an  advancing  crack 
causes  the  crack  surfaces  to  close  during  some  portion  of  loading  cycle.  He 
proposed  that  crack  growth  is  effected  by  the  portion  of  loading  cycle  over 
which  the  crack  remains  open.  This  concept  has  been  widely  advocated  in 
explaining  the  behavior  of  short  crack  acceleration  phenomenon.  The 
two-dimensional  results,  however,  indicate  that  large  amount  of  plasticity 
with  loss  of  constraint  rather  than  residual  or  crack  closure  effects,  leads 
to  the  breakdown  in  LEFM  conditions  near  the  crack  tip.  The  short  crack 
effects  grow  with  increasing  ratios  of  applied  stress  to  yield  stress,  and  for 
decreasing  crack  tip  constraint. 
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The  growth  of  the  surface  or  through  crack  is  three-dimensional,  where 
the  region  near  the  free  surface  is  governed  by  plane  stress  conditions,  and 
at  the  interior,  the  behavior  is  governed  by  plane  strain  conditions.  The 
small  flaw  acceleration  problem  may  be  due  to  the  influence  of 
three-dimensional  surface  flaw  plasticity  effects,  an  influence  that  is  not 
modeled  by  two-dimensional  approaches.  However,  the  task  of  modeling  the 
three-dimensional  propagating  crack  is  substantial.  The  only  numerical  study 
reported  to  date  is  the  limited  investigation  of  Chermahini  [3]  using  the 
finite  element  method.  Unlike  the  two-dimensional  boundary  element  modeling, 
where  the  presence  of  the  crack  is  taken  care  of  by  a  fundamental  solution 
and,  therefore,  does  not  require  the  modeling  of  crack,  the  three-dimensional 
modeling  requires  modeling  of  the  crack  surface.  The  elastic  problems  solved 
by  using  the  BEM  as  reported  here  substantiates  the  advantage  of  the  method 
for  this  class  of  problems.  However,  the  need  to  have  small  plastic  cells  to 
model  the  yielded  region  near  the  surface  in  an  elastoplastic  crack  problem 
imposes  a  very  high  numerical  constraint.  While  this  limitation  has  been 
overcome  to  a  certain  extent  in  this  study,  a  complete  investigation  requires 
further  research  effort.  Some  numerical  results  of  the  three-dimensional 
elastoplastic  propagating  crack  are  presented  at  the  end  of  this  report. 
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2.0  TWO-DIMENSIONAL  ELASTOPLASTIC  FRACTURE  MECHANICS  MODELING 


2 . 1  Problem  Statement 

Figure  1  illustrates  the  correlating  role  of  the  elastic  stress  intensity 
factor  (K)  for  fatigue  crack  growth  problems.  Three  different  fatigue 
specimens  (A,  B,  C)  were  used  to  obtain  baseline  crack  growth  data:  a  standard 
ASTM  modified  compact  tension  specimen,  plus  two  center-cracked  specimens  of 
different  crack  sizes  and  thicknesses.  Two  sets  of  surface  crack  growth  data, 
taken  from  replicas  of  the  surface  cracks  on  a  full-size  component,  are 
plotted  in  Figure  1  for  smooth  sections  (D)  and  for  notches  (E) . 

The  fact  that  the  data  are  correlated  for  a  wide  range  of  crack  size, 
crack  shape  and  local  stress  field  confirms  the  usual  applicability  of  linear 
elastic  fracture  mechanics  (LEFM)  for  many  problems.  Other  data,  e.g.,  [4-7] 
clearly  indicate  that  the  elastic  stress  intensity  factor  from  LEFM  does  not 
correlate  all  crack  data  for  problems  with  high  applied  stress,  relative  to 
the  material  yield  stress.  The  crack  sizes  for  these  data  are  quite  small 
compared  to  the  crack  sizes  used  for  laboratory  characterization  of  crack 
growth . 

Thus,  there  is  a  clear  indication  from  the  experimental  data,  that  LEFM 
mav  apply  over  a  wide  range  of  geometric  conditions  but  not  over  all 
geometries.  A  recent  review  by  Leis,  et  al.  [4]  attributes  the  so-called 
small  flaw  modeling  problem  to  a  lack  of  crack  tip  similitude.  Stated 
somewhat  more  broadly,  the  damage  process  in  the  critically  stressed  (or 
strained)  volume  of  material  ahead  of  the  crack  is  no  longer  independent  of 
scaling  below  a  certain  size,  which  will  be  seen  to  depend  on  the  applied 
stress  level.  Those  who  argue  for  the  loss  of  similitude  as  the  reason  for 
the  loss  of  correlation  in  the  small  flaw  regime  must  seek  a  new  measure  of 
damage  within  the  critical  material  volume,  a  measure  that  reduces  to  LEFM  for 
the  appropriate  class  of  conditions. 
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Figure  1.  Crack  Growth  Rate  Data  for  Surface  Cracked 
Coaponents  Coapared  to  Baseline  Data 
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The  presence  of  a  plastic  zone  ahead  of  all  fatigue  cracks  is 
generally  ignored  except  for  problems  of  overload  and  load  ratio  (R)  effects. 
Baseline  fatigue  crack  growth  (FCG)  data  is  usually  generated  at  low  load 
ratios  (R  -  0.1);  this  data  intrinsically  contains  any  effects  of  the  plastic 
zone.  Generally,  for  LEFM  applicability,  the  strain  field  surrounding  the 
plastic  zone  is  controlled  by  the  elastic  strain  field,  and  is  not  a  separable 
influence . 

The  plastic  zone  becomes  dominant  for  fatigue  crack  growth 
problems  where  load  history  effects  are  important  [8].  An  overload  is 
sometimes  seen  as  a  modification  of  the  residual  stresses  due  to  plasticity 
ahead  of  the  crack  which  changes  the  cyclic  stress  intensity  factor  range,  or 
a  change  in  the  range  of  loading  for  which  the  crack  tip  is  closed  (no  fatigue 
damage)  [9-11].  The  actual  solutions  reported  herein  will  demonstrate  the 
complete  duality  of  the  residual  stress  and  crack  closure  arguments. 

The  plastic  zone,  which  remains  along  the  crack  surfaces  after  the 
crack  tip  advances,  is  referred  to  as  the  plastic  wake.  Under  cyclic 
conditions,  the  tensile  stretch  of  the  material  ahead  of  the  crack  tip  exceeds 
the  amount  of  reversed  plastic  strain  as  the  load  is  released.  The  excess  of 
plastic  strain  acts  as  an  interference  at  the  original  crack  plane,  resulting 
in  crack  surface  closure,  even  at  positive  loading.  If  there  is  an  additional 
effect  on  the  crack  surface  of  interference  due  to  mechanical  locking  of 
irregular  crack  surfaces,  or  due  to  significant  oxidation,  then  the  crack  tip 
may  be  wedged  open  under  zero  load. 

The  actual  fatigue  process  in  the  near  crack  tip  region  is,  of 
course,  quite  complex.  Cyclic  plasticity  may  be  distributed  (wavy  slip)  or 
may  be  concentrated  (planar  slip).  The  structure  of  the  material  may  play  a 
key  role  in  the  crack  growth  kinetics  through  grain  size  or  grain  boundary 
resistance  to  crack  extension.  The  conditions  at  the  crack  tip  vicinity  in 
regard  to  these  phenomena  may  include  cyclic  strain  and  mean  stress  effects. 


Thus,  an  understanding  of  the  kinematics  of  the  deformation  and  the  state  of 
stress  are  likely  to  be  es  .ential  to  an  understanding  of  the  large  crack 
versus  short  crack  growth  rate  data. 


2 . 2  Crack  Tip  Singular  Fields 


2.2.1  Review  of  BIE  for  Elastoolastlc  Cracks 

The  elastoplastic  boundary -integral  equation  (BIE)  method  is  used 
for  the  current  study.  The  formulation  of  the  BIE  uses  the  elastic  Green  s 
function  to  represent  a  planar  body  with  a  straight  line,  stress  free  crack 
[12] .  The  plasticity  is  a  classical  von  Mises  plasticity  formulation,  as 
recently  reported  [13,14].  The  resulting  elastoplastic  BIE  formulation  for 
the  fracture  mechanics  problem  has  yielded  some  new  and  important  insight  into 
the  physical  problem  and  will  be  briefly  reviewed. 


Figure  2  shows  a  finite  plate  with  a  central  or  edge  crack  which 
is  taken  to  be  stress  free.  The  plastic  strains  are  assumed  to  be  contained 
within  the  discrete  element  region  around  the  end  of  the  crack;  the 
elastoplastic  strains  are  taken  to  be  piecewise  constant  within  each  of  the 
domain  elements.  Using  the  fundamental  solution  for  a  point  load  in  an 
infinite  elastic  plate  containing  a  single  straight  crack  of  length  2a,  the 
Somigliana  identity  for  the  interior  strain  increment  fjj  for  the 
elastoplastic  problem  is  given  from  [13]: 
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In  Eq.  (1),  Djjjc  and  Sjjk  are  the  symmetric  parts  of  the  gradients  of  the 
boundary  tractions  (on  the  uncracked  surface  S)  and  boundary  displacements  due 
to  a  unit  load  at  an  arbitrary  interior  point.  The  tensor  G^ij  are  generally 
the  symmetric  part  of  the  gradient  of  the  interior  stress  tensor  due  to  the 
elastic  load.*  This  tensor  kernel  is  highly  singular  and  has  been 

*  This  is  true  for  three-dimensional  elastoplasticity  and  for  the  plane 
stress  kernel.  The  plane  strain  kernel  is  of  a  different  form  [13]. 


Figure  2.  Domain  Discretization  for 
Long  and  Short  Cracks — Stationary  Crack 


differentiated  at  the  singular  point  boundary  which  is  excluded  from  the 
integral  in  a  principal  value  sense.  As  a  result  of  the  differentiation  at 
the  boundary  of  the  exclusion  region  surrounding  the  singularity,  one  obtains 
Eiklj ■  an  elastic  constant  tensor  which  differs  for  plane  strain  and  plane 
stress.  The  boundary  data  is  given  as  the  increment  in  boundary  traction  t. 
and  boundary  displacement  u  ;  the  increment  in  plastic  strains  (inside  plastic 
domain  A)  is  given  by 

Equation  (1)  should  be  seen  as  the  complete  statement  regarding 
the  elastoplastic  field  problem.  That  is,  the  total  strain  field  is 
mathematically  complete.  The  stresses  obtained  from  (1)  can  be  shown  to  fully 
satisfy  equilibrium  at  all  interior  points.  The  requirement  on  the  plastic 
strain  increments  is  only  that  they  be  no  more  singular  than  a  1/R 
singularity,  where  R  is  the  distance  from  the  crack  tip,  (true  for  perfect 
plasticity  at  crack  tips),  as  discussed  in  [13].  The  plastic  strains  are 
allowed  to  be  discontinuous  and  permit  a  slip  line  interpretation  for  fully 
developed  plasticity  problems. 

2.2.2  Strain  Intensity  Factor 

The  total  strain  field  in  Eq.  (1)  contains  both  explicit  and 
implicit  singular  behavior  near  the  crack  tip.  The  details  of  the  singular 
behavior  are  discussed  in  detail  in  [13]  and  are  now  briefly  outlined.  The 
boundary  terms  in  Eq.  (1)  explicitly  depend  on  the  inverse  of  the  square-root 
of  the  crack  tip  distance  near  the  crack  tip.  The  influence  of  plastic  strain 
is  discussed  in  [11,12]  for  both  singular  and  non- singular  plastic  strains. 

When  the  plastic  strains  are  non-singular,  then  the  total  strain 
contribution  from  the  integral  in  Eq.  (1)  is  also  inversely  proportional  to 
the  square-root  of  the  crack  tip  distance.  For  such  problems  the  total  strain 


contains  the  same  singular  behavior  as  for  the  LEFM  problem.  A  stress 
intensity  factor  result,  based  on  the  strain  singularity  is  given  from  [10] 
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In  Eq.  (2)  the  mixed  mode  (Mode  I  and  Mode  II)  stress  intensity  factors  are 
given  as  simple  vector  products  of  the  boundary  data  with  the  boundary  vector 
potentials,  and  L^,  and  the  inner  product  of  the  total  plastic  strains  with 
the  volume  potential  tensor  M^j ,  taken  over  the  domain  of  plastic  strains. 

For  the  case  of  singular  plastic  strains,  Eq.  (1)  reduces  to  a 
homogeneous  integral  equation  containing  only  the  plastic  singularity.  The 
order  of  the  singularity  is  shown  in  [10]  to  be  implicitly  determined,  based 
on  the  numerical  results  for  the  strain  hardening  characteristic  of  the 
material.  The  limit  on  the  order  of  the  singularity  for  convergence  of  the 
volume  integral  in  Eq.  (1)  has  been  shown  to  correspond  to  the  plastic  strain 
singularity  for  perfect  plasticity  (1/R). 


Equation  (2)  provides  a  unique  means  for  investigating  the  plastic 
zone  of  a  stationary  or  propagating  fatigue  crack.  At  any  state  of  plastic 
strains,  Eq.  (2)  may  be  used  to  compute  a  strain  intensity  factor  at  the  crack 
tip.  The  plastic  strain  distribution  for  a  given  crack  size  affects  the 
residual  stresses  ahead  of  the  crack.  This  residual  stress  effect  may  be 
calculated  from  Eq.  (2)  by  increasing  the  crack  length  parameter  "a"  to  a  +  Aa, 
and  then  by  allowing  an  elastic  iteration  of  the  numerical  solution;  the 
elastic  iteration  is  required  to  satisfy  the  traction- free  boundary  conditions 
on  the  new  crack  surface.  The  procedure  is  entirely  equivalent  to  calculating 
the  residual  stress  effect  for  any  elastic  crack  propagating  through  a  plastic 
strain  distribution  due  to  prior  load  history  (see  [15]).  In  the  case  of  the 
plastic  strains  due  to  the  crack  itself,  Eq.  (2)  has  been  analytically  shown 
[13]  to  be  nonconvergent  at  da  ”  0!  numerical  data  in  [15]  confirmed  this.  The 


interpretation  at  Aa>o,  Eq.  (2)  provides  a  direct  means  for  calculating  the 
effect  of  all  prior  plastic  strains  on  the  crack  tip  stress  intensity  factor. 
It  should  be  noted  that  while  Eq.  (2)  is  reported  in  terms  of  stresses,  it  is 
directly  computed  from  the  strain  field  limit  approaching  the  crack  tip  times 
the  material  stiffness  matrix.  Thus,  the  result  is  in  reality  a  strain 
intensity  factor.  The  stresses  have  no  singular  interpretation  for 
elastoplastic  fracture  mechanics  analysis. 

The  direct  application  of  this  mathematical  condition  is  to 
calculate  retardation  effects  and  crack  closure  effects  for  fatigue  crack 
overload  conditions.  The  details  of  the  overload  retardation  and  closure 
calculations,  as  well  as  an  exploration  of  the  implications  for  fatigue  crack 
growth  are  reviewed  in  the  following  sections  of  this  text. 

2.2.3  Crack  Opening  Displacement 

The  Somigliana  identity  for  the  interior  displacements  in 
elastoplastic  fracture  mechanics  problem  is  given  as: 

ui  3  *  T1jujdS  +  £  U1jVS  +  <^>r1Jk€jkdA  <3> 

where  Tjj  and  Ujj  are  the  boundary  tractions  and  displacements  on  the 
uncracked  surface  S  for  point  loads  applied  in  each  of  the  xj  directions.  The 
tensor  I^j|<  is  generally  the  interior  stress  state  for  the  elastic  cracked  body 
due  to  the  point  loads.  The  displacement  in  Eq.  (3)  may  be  evaluated  on  the 
boundary  of  the  crack,  without  major  algorithmic  problems.  It  is  easily 
shown,  following  the  same  procedures  as  in  [13],  that  the  crack  opening  field 
near  the  crack  tip  contains  a  square- root -R  behavior  due  to  the  boundary  data, 
and  a  term  due  to  the  (singular)  plastic  strains. 

The  crack  opening  displacement  result  exactly  parallels  the  one 
above  for  the  total  strains.  Thus,  the  effect  of  the  plastic  strains  on  crack, 
opening  displacement  is  directly  linked  to  the  effect  of  plastic  strains  on 
the  strains  ahead  of  the  crack  tip.  Thus,  if  the  plastic  strains  are  bounded, 


or  Che  elastic  crack  tip  is  extended  to  a+Ad ,  the  total  crack  opening  field 
near  the  crack  tip  is  dependent  on  the  square-root  of  the  crack  tip  distance, 
R.  For  the  elastoplastic  fracture  mechanics  solution,  the  near  crack  tip 
behavior  will  depend  on  the  implicit  crack  tip  singularity  of  plastic  strains. 
The  numerical  results  will  expand  on  this  more  fully. 

2.2.4  Plastic  Slip  Line  .Solutions 

The  last  area  of  potential  exploitation  of  the  elastoplastic  BIE 
formulation  for  the  fracture  mechanics  problem  is  its  use  in  modeling  fully 
developed  plastic  fields.  As  discussed  in  [16],  the  elastoplastic  field  is 
theoretically  hyperbolic,  such  that  the  plastic  strains  may  be  discontinuous 
along  characteristic  lines  whose  directions  depend  on  the  solution  geometry. 

The  BIE  formulation  in  Eq.  (1)  admits  such  discontinuous  plastic 
strains  without  difficulty.  The  derivation  of  Eq.  (1)  requires  only  that  the 
plastic  strains  be  integrable  in  the  volumetric  sense;  no  derivatives  of  the 
plastic  strains  are  required.  The  boundary  terms  in  Eq.  (1)  are  the  result  of 
an  integration  of  the  displacement  gradients  through  the  divergence  theorem. 
The  application  of  the  divergence  theorem  only  requires  that  the  displacements 
be  continuous  and  differentiable,  a  condition  satisfied  by  compatibility. 

An  example  of  the  application  of  Eq.  (1)  and,  more  specifically, 

Eq.  (2)  to  a  problem  of  a  discontinuous  plastic  field  was  given  in  [12]  for 
the  problem  of  a  modeled  weld  bead  with  a  narrow  zone  of  high  residual  plastic 
strain.  The  application  was  to  the  computation  of  crack  tip  stress  intensity 
factors  for  cracks  within  the  plastic  field,  and  for  crack  tips  extending 
beyond  the  plastic  field.  The  solution  was  essentially  exact. 

The  implications  of  this  result  are  not  yet  fully  explored. 

However,  it  is  clear  that  the  BIE  formulation  for  this  class  of  problems 
admits  the  possibility  of  modeling  slip  line  plastic  strain  solutions  without 
the  need  for  special  domain  elements  or  special  element  refinement  near  the 


zone  of  discontinuity.  This  contrasts  with  the  finite  element  method  which 
requires  the  use  of  many  small  elements  along  such  lines  of  discontinuity  for 
any  degree  of  accuracy. 


The  material  is  assumed  to  be  elastic-perfectly  plastic.  The 
modulus  is  taken  to  be  30*10)6  pSi  while  the  yield  stress  is  taken  to  be  50 
ksi.  These  are  arbitrary  constants  and  none  of  the  conclusions  specifically 
depends  on  these  values.  The  stationary  load  case  refers  to  a  load/unload 
cycle  at  the  original  crack  length.  Cyclic  crack  modeling  refers  to  a 
load/unload  cycle,  a  crack  extension  by  some  increment,  and  then  another 
load/unload  sequence,  etc.  A  validation  model  for  the  load  application  phase 
using  the  elastoplastic  BIE  code  is  reported  in  [14],  and  is  repeated  in 
Figure  3 .  The  results  were  compared  in  Figure  3  to  numerical  FEM  results  from 
ADINA  to  establish  the  validity  of  the  results. 

Two  crack  lengths  are  used  in  the  current  study,  based  on  a  desire 
to  have  small-scale  and  large-scale  yielding.  The  crack  lengths  are  taken  to 
be  10.0"  and  0.1"  with  applied  stress  intensity  factors  (elastic)  of  20  ksi. 
/in.  The  plate  is  taken  to  be  1000"  wide  and  long  to  minimize  the  effects  of 
finite  geometry.  Based  on  a  Dugdale  model,  the  plastic  zone  size  for  plane 
stress  is  0.062"  for  both  crack  sizes.  The  plane  strain  value  might  be 
estimated  to  be  one- third  of  that  size,  or  0.021". 

Thus,  in  the  case  of  the  large  crack,  the  normalized  size  of  the 
plastic  zone  (rp/a)  is  0.0063,  while  for  the  short  crack  case  the  value  is 
0.63.  Clearly,  the  two  problems  do  not  possess  similitude  of  the  plastic 
strain  fields  and  may  serve  as  good  comparisons  for  large  crack  and  short 
crack  results. 
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Figure  2  shows  the  domain  discretization  currently  used  for 
analysis  of  the  stationary  (non- propagating)  crack.  The  mesh  is  much  cruder 
than  used  in  the  study  reported  previously  [13]  but  produces  essentially  the 
same  plastic  strain  results.  The  loading  was  applied  incrementally  to  an 
equivalent  elastic  stress  intensity  factor  of  20  ksi./in  and  thereby  removed 
incrementally.  The  plastic  strains  at  the  crack  tip  are  thereby  reversed,  but 
no  recycling  of  the  loading  is  performed. 

In  order  to  define  the  size  of  the  plastic  zone  in  an  unambiguous 
manner,  use  is  made  of  the  concept  of  elastic  crack  extension  through  the 
plastic  strains  for  the  stationary  crack  analysis.  For  this  analysis,  the 
stress  (strain)  intensity  factor  is  computed  using  Eq.  (2)  for  various  crack 
extensions;  the  plastic  zone  taken  from  that  for  the  the  original  crack  size 
and  the  above  loading  conditions  and  material  properties . 

Figures  4  and  5  show  the  resulting  values  of  elastic  stress 
intensity  factors.  As  discussed  in  [14],  the  computation  for  (2)  treats  the 
prior  load  plasticity  as  a  residual  strain  distribution  and  provides  the 
effect  of  the  residual  strains  on  the  stress  intensity  factor  for  various 
crack  lengths.  The  elastoplastic  solution  is  first  obtained  so  that  a  plastic 
strain  distribution  is  known.  Then  the  BIE  code  is  exercised  in  a  strictly 
elastic  sense.  The  crack  size  is  incremented  by"  4a"  into  and  through  the 
plastic  zone  of  the  previously  analyzed  elastoplastic  condition.  Each  elastic 
calculation  from  Eq.  (2)  is  then  a  computation  of  any  elastic  stress  intensity 
factor,  including  the  effect  of  the  prior  elastoplastic  zone.  The  resulting 
stress  intensity  factor,  while  not  strictly  physical,  directly  computes  the 
effect  of  the  previous  plasticity  on  the  crack  tip  strain  field  for  the  larger 
size  crack.  Thus,  this  is  an  unambiguous  estimate  of  the  residual  plastic 
zone  effect  for  propagating  fatigue  cracks.  The  closest  parallel  to  previous 
crack  modeling  strategies  is  the  estimate  of  overload  plasticity  effects  on 
the  elastic  stress  intensity  factor  for  the  fatigue  load  condition  [8]. 
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Figure  5.  Residual  Strain  Effect  on  Crack  Tip  Stress  Intensity  Factor 
for  Various  Crack  Extensions '-Short  Crack  at  No  Load 


Table  1  compares  the  estimated  plastic  zone  size  from  these 
results  to  the  Dugdale  (plane  stress)  and  plane  strain  plastic  zone  sizes. 

Table  1 

Comparison  of  Plastic  Zone  Sizes 


Long  Crack  rp(est. )  rp(cal. ) 

Plane  Stress  0.063  0.060 

Plane  Strain  0.021  0.025 

Short  Crack 

Plane  Stress  0.063  0.110 

Plane  Strain  0.021  0.080 


As  can  be  seen,  the  BIE  results  compare  quite  favorably  to  the  estimates  for 
the  long  crack  case  but  not  for  the  short  crack  case.  The  effect  of  the 
plastic  strains  for  the  short  crack  case  is  substantially  greater  than  for  the 
long  crack  case;  the  difference  between  the  plane  strain  and  plane  stress 
plastic  zones  is  quite  diminished.  These  results  seem  quite  consistent  with 
the  expectation. 

2.3.2  Ctask  Shensi, <?n  Modeling 

Figure  6  shows  the  domain  discretization  for  the  plastic  strains 
used  for  a  study  of  the  crack  propagating  as  in  fatigue.  That  is,  the  crack 
is  incrementally  loaded  to  20  ksi./in  and  unloaded  at  each  crack  size.  The 
yield  stress  in  compression  is  taken  to  be  the  same  as  in  tension  (no 
Bauschinger  effect) .  The  original  crack  tip  is  as  shown  for  the  length  "a" . 

In  the  current  study,  crack  tip  plastic  strains  are  allowed  to  reverse  before 
the  crack  is  extended.  This  strain  reversal  accounts  for  what  others  model  as 
the  Interference  of  the  plastic  wake  for  a  propagating  crack  [4,6]. 

The  surfaces  of  the  crack  are  not  constrained  against 
inter-penetration  upon  further  unloading  to  zero  load.  By  reversing  the 
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plastic  strains  near  the  crack  tip  during  unloading  and  before  the  crack  has 
been  extended,  the  primary  effect  of  crack  surface  interaction  on  the  plastic 
strains  in  the  wake  is  computed.  Somewhat  higher  cyclic  plastic  strains  are 
predicted  very  near  the  crack  tip,  due  to  the  additional  plastic  strain 
occurring  at  loads  lower  than  the  load  at  which  closure  occurs.  The  bulk  of 
the  closure  effect  beyond  reversing  the  plastic  strains  for  the  material  just 
ahead  of  the  crack  tip  prior  to  extension  is  the  elastic  problem  of  loads  on 
the  crack  surface  to  maintain  zero  closure.  This  nonlinear  contact  stress 
problem  is  fully  reversible  upon  load  application  and  does  not  affect  the 
plastic  results. 

The  crack  length  is  increased  after  unloading  so  the  crack  tip 
will  be  at  the  next  domain  element  boundary  (0.0125  units  long  in  the  domain 
mesh  used) .  Extension  of  the  crack  is  achieved  in  the  elastoplastic  BIE  code 
by  assigning  a  new  length  to  "a"  and  then  by  a  single  iterative  cycle  to 
achieve  elastic  load  redistribution.  Thus,  the  traction  free  condition  on  the 
crack  surface  is  exactly  satisfied  without  any  mesh  changes.  The  load/unload 
cycle  is  repeated,  followed  by  crack  extension,  etc.  until  the  crack 
plasticity  results  stabilize.  The  stabilization  of  the  cyclic  plastic  zone 
took  about  2  to  3  crack  extensions.  Stability  of  the  numerical  results  was 
greater  for  the  plane  strain  than  for  the  plane  stress  condition. 

Figures  7  and  8  plot  the  results  obtained  by  calculation  of  the 
stress  intensity  factor  from  Eq.  (2)  for  various  elastic  crack  tip  locations 
displaced  from  the  tip  of  the  analyzed  crack  size  ("a  +  J  a").  As  in  the  case 
of  the  stationary  crack  reported  above,  the  computation  of  stress  intensity 
factor  is  a  direct  calculation  of  the  effect  of  prior  plasticity  on  a 
propagating  fatigue  crack. 

In  all  problems  only  three  extensions  of  the  crack  showed  that  the 
plastic  strains  had  converged  acceptable  to  a  steady-state  condition, 
particularly  for  the  long  crack  solution.  In  fact,  the  cyclic  results  did  not 
differ  much  from  the  monotonic  results.  The  plane  stress  results  for  the 
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Figure  8.  Residual  Strain  Effect  on  Crack  Tip  Stress  Intensity  Factor 
for  Propagating  Cracks  at  No  Load — Short  Crack 


short  crack  showed  the  weakest  convergence,  both  In  terms  of  the  crack 
extension  comparisons  but  also  in  terms  of  the  numerical  algorithm.  The  short 
crack,  plane  stress  case  has  the  weakest  constraint  on  the  plastic  strains, 
and  this  would  seem  a  reasonable  source  of  the  difficulty. 

It  can  be  seen  that  the  computed  plastic  stress  intensity  factors 
are  somewhat  oscillatory  for  short  crack  extensions.  The  oscillations  occur 
when  the  crack  extension  is  inside  the  smallest  crack- tip  domain  element. 
However,  by  ignoring  the  plastic  strain  (Eq.  (2))  within  the  smallest  domain 
element,  the  elastic  stress  intensity  factor  is  made  regular  as  seen  in  Figure 
9.  The  Figure  shows  increasingly  negative  values  of  strain  intensity  factor 
from  Eq.  (2)  as  the  smallest  domain  element  is  reduced  in  size.  This 
divergence  in  K  occurs  as  expected  due  to  the  higher  order  plastic  singularity 
in  crack  tip  strain.  Vhat  was  not  expected  was  the  apparently  small  (perhaps 
indefinitely  small)  zone  over  which  the  divergence  affects  the  calculated  K  in 
Eq.  (2).  This  zone  was  seen  to  reduce  as  the  domain  element  as  the  crack  tip 
was  reduced.  This  implies  a  discontinuity  in  the  behavior  of  Eq.  (2)  for 
small  da.  At  da-o,  the  result  is  divergent.  For  da> o,  the  result  seems  to 
converge  to  finite  values  of  K  for  small  da. 

The  cross -data  points  in  Figure  9  correspond  to  values  from  Eq. 

(2)  obtained  by  setting  the  plastic  strain  in  the  smallest  element  to  zero. 

It  is  seen  in  Figure  9  that  the  effect  of  the  plastic  strain  in  the  first  and 
smallest  domain  element  zone  at  the  crack  tip  is  limited  to  da  values  only 
within  that  element. 

2.3.3  Closure  vs.  Retardation  Results 

The  crack  tip  opening  was  calculated  from  Eq.  (3)  for  both 
stationary  (monotonic  loading)  and  propagating  cracks.  The  algorithm  is  based 
on  monitoring  the  displacement  at  a  pre- defined  point  on  the  crack  surfaces  at 
a  distance  of  .00001"  behind  the  crack  tip.  This  is  an  arbitrary  location  but 


reflects  the  experimental  capability  [11].  Similitude  is  not  maintained  as 
0.0001/a  equals  10 ‘6  for  the  long  crack  case  and  equals  10 for  the  short 
crack  case. 

Loading  for  crack  closure  monitoring  is  applied  elastically,  as 
the  plastic  effect  cannot  resume  until  after  the  crack  tip  is  open.  The  load 
at  which  the  displacement  at  this  monitoring  location  goes  from  negative 
(overleaving  of  the  crack)  to  positive  is  defined  as  the  crack  tip  opening 
load.  Interpolation  of  the  numerical  results  was  used  to  select  this  point. 
The  load  for  opening  is  then  converted  into  an  equivalent  elastic  stress 
intensity  factor  corresponding  to  crack  opening. 

The  results  for  the  propagating  crack  are  shown  in  Figures  10  and 
11.  The  methodology  used  for  the  residual  solution  is  obtained  for  a  given 
crack  size,  then  the  BIE  code  elastically  computes  the  effect  of  the  plastic 
strains  on  the  crack  opening  load  for  various  elastic  extensions  of  the  crack. 
This  corresponds  then  to  calculating  the  effect  of  the  plastic  wake  on  crack 
closure.  It  should  be  recalled  that  the  elastoplastic  loading  and  unloading 
are  done  for  a  given  size  crack  so  that  the  effects  of  interference  of  the 
plastic  wake  from  both  sides  of  the  crack  is  taken  fully  into  account  (see 
[7]). 

The  numerical  data  in  Figures  7-8  and  10-11  show  a  very  important 
piece  of  insight  into  the  fatigue  crack  growth  problem.  Namely,  the  effect  of 
the  plastic  wake  on  the  apparent  stress  intensity  factor  for  crack  opening  is 
identical  (with  a  sign  change)  to  the  effect  of  the  residual  stress  ahead  of 
the  crack  on  retardation  of  the  stress  intensity  factor.  Apparently,  Kanninen 
and  his  covorkers  reported  the  same  finding  for  the  super -dislocation 
plasticity  model  discussed  in  Chapter  8  of  his  book  [6],  The  theoretical 
reason  for  this  is  clearly  seen  from  Eqs.  (1)  and  (3),  which  relate  to  crack 
tip  strain  field  and  crack  opening  displacements.  The  strain  field  Eq.  (1)  is 
computed  by  differentiation  of  the  displacement  field  Eq.  (3)  resulting  in  a 
mathematical  equivalence  for  any  elastoplastic  field. 
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Figure  11.  Residual  Strain  Effect  on  Crack  Tip  Opening  Load — Short  Crack 


The  duality  of  residual  stress  and  crack  closure  effects  was  also 
found  for  the  case  of  the  stationary  crack  subjected  to  an  elastoplastic 
load/unload  solution.  The  application  of  this  result  would  be  for  overload 
modeling,  with  the  conclusion  that  accurate  modeling  of  retardation  and 
closure  effects  of  the  overload  would  produce  identical  effects  on  the  crack 
tip  cyclic  response.  This  is  true  during  fatigue  crack  propagation  with 
plastic  wake  effects,  as  crack  tip  opening  occurs  as  the  last  event  after  the 
rest  of  the  crack  has  been  opened  during  loading  [6,7], 

2.3.4  Crack  Opening  Displacement 

Figure  12  plots  normalized  crack  opening  displacements  for  the 
long  crack  under  plane  strain  conditions  at  the  crack  length  of  a-10.0.  The 
displacements  are  normalized  by  the  elastic  crack  opening  displacement 
function: 

»  *  B7T7  (2«r-r2)1/2  (*> 
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where  r  is  the  crack  tip  distance  and  H*E,  E/(l-v)  for  plane  stress  and  plane 
strain.  The  solutions  are  plotted  during  loading  unloading  and  show  the  crack 
opening  during  loading  at  K  -  5  ksi./in,  which  is  elastic  due  to  the  finite 
size  of  the  smallest  domain  element  for  plastic  strain  modeling,  at  K  -  10 
ksi./in  and  at  20  ksi./in.  The  exact  solution  is  seen  to  fit  the  elastic 
behavior  of  the  crack  tip  solution  out  to  a  value  of  r/a  -  0.01.  At  higher 
stress  intensity  factors,  the  crack  opening  displacement  is  seen  to  be 
increasingly  compliant,  resulting  in  higher  displacements  over  a  length  from 
the  crack  tip  somewhat  larger  than  the  plastic  zone  size  of  the  problem. 

During  unloading,  only  a  portion  of  the  plastic  stretch  is  reversed,  thus,  the 
ratio  of  computed  to  elastic  crack  opening  displacement  becomes  progressively 
higher.  The  same  behavior  was  observed  for  both  plane  strain  short  crack  and 
plane  stress  long  and  short  crack  cases. 


Figure  13  plots  the  same  results  for  the  case  of  the  crack  which  has  been 
extended  by  0.0025a.  The  plastic  wake  effect  has  the  result  of  reducing  the 
crack  opening  displacements  to  the  elastic  values  down  to  r/a  <  0.001.  Given 
a  reasonable  hypothesis  that  cyclic  crack  extension  is  associated  with  crack 
tip  cyclic  deformation,  then  the  cyclic  crack  opening  displacement  values  for 
the  long  crack  plane  strain  results  indicate  that  the  LEFM  model  should  be 
applicable.  Of  course,  this  is  just  as  expected.  The  role  of  the  elastic 
strain  field  for  the  small-scale  yielding  problem  is  to  force  the  plastic  zone 
to  conform  to  the  deformation  of  the  elastic  singularity,  near  the  crack  tip. 
The  conformance  is  clearly  achieved  during  the  unloading  when  reversed 
plasticity  is  incurred. 

Figure  14  plots  the  same  type  of  COD  data  for  the  short  crack 
problem  where  the  crack  is  extended  by  0.25a.  The  effect  of  plasticity  on 
crack  opening  compliance  is  now  seen  by  the  reduction  of  the  COD  below  the 
elastic  value,  thus,  achieving  closure  at  positive  loading.  Therefore,  the 
crack  opening  displacement  values  for  short  crack  results  indicate  that  the 
LEFM  model  is  not  applicable.  The  data  also  shows  large  relative  crack 
opening  near  the  crack  tip,  indicative  of  blunting. 

The  plane  stress  COD  solution  for  both  long  and  short  cracks  are 
plotted  in  Figures  15  and  16,  respectively.  The  plane  stress  results  show  an 
enduring  plastic  wake  for  both  long  and  short  cracks;  the  difference  between 
the  plane  strain  and  plane  stress  results  is  decreased.  The  increased  plastic 
wake  size  for  plane  stress  is  due  to  diminished  constraint.  It  should  be 
noted  that,  realistically,  the  plane  stress  condition  occurs  only  at  very  thin 
plates  or  at  a  very  small  layer  at  the  boundary.  Similitude  for  long  crack 
fracture  mechanics,  therefore,  implies  that  plane  strain  elastoplastic 
behavior  must  be  occurring  for  real,  three-dimensional  cracks. 
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The  two  model  problems  considered  in  this  study  are  for  equally 
applied  elastic  stress  intensity  factors.  In  the  short  crack  case,  this 
results  in  a  ratio  of  applied  stress  to  yield  stress  of  slightly  more  than 
71%.  At  this  load  level,  the  elastoplastic  plane  strain,  long  crack  results 
predict  that  the  crack  responds  in  accordance  with  LEFM. 

The  short  crack  results  show  that  the  effect  of  plastic  wake  on 
COD  at  a  region  near  the  crack  tip  is  predicted  to  be  significant  for  the 
plane  strain  and  plane  stress  cases.  The  plastic  wake  for  short  cracks  is 
seen  to  be  the  result  of  high,  unreversed  crack  tip  stretch  plasticity.  The 
fact  that  the  crack  tip  stretch  within  the  scale  of  the  cyclic  crack  tip  zone 
is  greater  than  the  LEFM  prediction  may  provide  the  mechanics  for  increased 
short  crack  fatigue  crack  growth.  The  increased  crack  tip  stretch  shows  up 
physically  as  a  plastic  wake. 

2.4  Closure  and  Growth  of  Cracks 

One  of  the  important  "spects  of  the  fatigue  crack  analysis  is  to 
correlate  and  predict  the  rate  of  growth  of  cracks  for  given  combinations  of 
geometry,  material  properties,  and  load  history.  Many  theories  have  been 
proposed  to  predict  the  crack  growth  rate  of  cracks  that  fall  under  the  small 
scale  yield  categories.  Most  of  the  current  models  are  based  on  the  concept 
that  average  growth  per  loading  cycle  can  be  related  as  a  power  low  function 
of  the  stress  intensity  factor,  K  [17,18],  However,  under  large  scale 
yielding  conditions,  accurate  prediction  of  crack  growth  is  not  possible  even 
with  the  use  of  modified  forms  of  stress  intensity  factors.  Many  alternative 
parameters  such  as  strain  based  stress  intensity  factors  (e.g.,  [19,20]), 
crack  opening  displacements  (e.g.,  [21-24]),  equivalent  stress  intensity 
factors  (e.g.,  [25-27]),  and  J-Integrals  (e.g.,  [28,29])  have  been  suggested 
for  the  investigation  of  fatigue  crack  growth  where  LEFM  is  not  applicable. 
However,  since  first  noted  by  Elber,  crack  closure  concept  is  considered  as  a 
very  important  criterion  in  explaining  the  propagation  of  the  crack.  He  used 


an  effective  stress  intensity  factor,  corresponding  to  the  portion  of  loading 
over  which  the  crack  remains  open,  as  a  correlating  parameter.  In  this 
section,  we  further  study  the  effect  of  crack  closure  on  crack  growth.  A 
detailed  study  of  crack  closure  under  plane  stress  conditions  using  finite 
element  method  has  been  reported  recently  by  McClung  [30]. 

McClung's  investigation  focussed  on  crack  closure  and  factors  that 
control  its  magnitude  for  constant  amplitude  of  stress  level.  He  noted  the 
pronounced  decrease  in  normalized  opening  stress  with  increasing  maximum 
stress  that  is  largely  ignored  in  many  experimental  research.  Also,  he 
confirmed  an  earlier  finding  by  Iyyer  and  Dowling  [31,32]  that  the  ratio,  U, 
of  the  effective  stress  over  which  the  crack  remains  open  to  the  total  applied 
stress  depends  on  the  maximum  stress  amplitude. 

Five  crack  lengths  ranging  from  0.07  in.  to  10  in.  were  selected  to 
examine  the  dependence  of  U  on  the  maximum  stress  amplitude  in  the  current 
study.  In  addition  to  the  constant  stress  amplitude  case,  S(max) , 
investigated  by  McClung,  crack  closure  due  to  constant  stress  intensity  factor 
amplitude  case,  K(max) ,  under  plane  stress  as  well  as  plane  strain  conditions 
were  studied.  For  K(max)  case,  the  value  of  maximum  stress  was  reduced 
progressively  with  crack  advance.  The  values  of  the  ratio  of  maximum  stress 
amplitude  to  the  yield  stress  at  the  initial  crack  length  for  various  cases 
considered  are  given  below. 
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As  explained  in  Section  2.3,  the  crack  was  extended  by  arbitrary 
increments  of  0.0125  in.  at  zero  load  and  incrementally  loaded  to  the  maximum 
value  and  unloaded,  subsequently.  In  the  first  study,  the  value  of  K(max)  was 
kept  constant.  As  noted  previously  in  Figures  14-16,  the  crack  remained 
closed  during  a  portion  of  the  loading  cycle.  The  crack  closing  stresses  were 
computed  by  monitoring  the  maximum  load  at  which  the  crack  surfaces  initiate 
contact  during  the  unloading  cycle.  It  was  generally  observed  that  the 
initial  contact  occurred  at  or  near  the  location  of  last  crack  tip  before  the 
current  extension.  The  crack  closing  stress  normalized  with  respect  to  the 
maximum  stress  amplitude  is  plotted  during  crack  extension  in  Figure  17  for 
the  plane  stress  case  and  in  Figure  18  for  the  plane  strain  case.  The  results 
indicate  that  after  an  initial  increase,  the  crack  closing  loads  stabilize,, 
generally  after  about  5  increments. 

To  see  the  dependence  of  the  crack  opening  stress  on  the  maximum 
stress  amplitude,  the  stable  crack  closing  stresses  normalized  with  respect  to 
the  maximum  stress  for  different  maximum  stress  amplitudes  are  plotted  in 
Figure  19  for  the  plane  stress  case  and  in  Figure  20  for  the  plane  strain 
case.  The  results  indicate  that  the  normalized  crack  closing  stress  increases 
with  increase  in  maximum  stress  amplitude.  In  other  words,  the  normalized 
crack  closing  stress  increases  with  decreasing  crack  length,  going  from  long 
crack  to  short  crack.  These  findings  are  inconsistent  with  the  finite  element 
method  crack  opening  stress  results  reported  for  the  plane  stress  case  using 
constant  amplitude  loading  [30,33]. 

To  examine  the  difference  in  the  crack  opening  and  closing 
stresses,  the  crack  was  reloaded  at  the  final  crack  length  without  crack 
extension.  A  comparison  of  the  normalized  crack  opening  and  closing  stresses 
shown  in  Figure  21  indicates  that  whereas  the  crack  opening  stresses  are 
higher  than  the  closing  stresses  for  the  plane  stress  case,  the  values  are 
essentially  identical  for  the  plane  strain  case.  However,  the  crack  opening 
stresses  still  show  increase  with  increase  in  applied  maximum  load. 
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Figure  18.  Crack  Closing  Stresses  Through  Crack  Extension 
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Figure  20.  Crack  Closing  Stresses  at  Stable  State 
Plane  Strain,  Constant  K  (sax) 
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The  same  investigation  was  conducted  at  constant  stress  loading. 
The  normalized  crack  closing  stresses  for  various  crack  extensions  are  plotted 
in  Figure  22  for  the  plane  stress  case  and  in  Figure  23  for  the  plane  strain 
case.  The  results  seun  to  stabilize  after  about  10  crack  increments,  however, 
the  stabilization  is  weaker  for  this  loading  than  the  constant  K  loading 
reported  earlier,  especially  for  short  crack  cases.  The  normalized  crack 
closing  stresses  for  various  stress  amplitudes  for  the  plane  stress  case  shown 
in  Figure  24  indicate  that  the  closing  stresses  decrease  with  increasing 
applied  maximum  loads,  which  is  consistent  with  the  results  reported  by 
McClung  [30]  and  Ogura,  et  al.,  [33]  for  high  stress  amplitudes.  The 
plane  strain  results  in  Figure  25  also  indicate  decrease  in  closing  stresses 
at  high  stress  amplitudes. 

It  should  be  noted  that  there  are  discrepancies  in  the  modeling 
approaches  used  by  the  finite  element  method  [30,33]  and  the  current  BEM 
approach.  In  the  finite  element  method  the  crack  extension  was  simulated  by  a 
node  release  scheme  at  the  maximum  load  level  and  the  crack  opening  stresses 
were  computed  during  the  subsequent  loading  cycle.  In  the  current  approach 
the  crack  was  extended  to  zero  load  and  the  closing  levels  were  monitored 
during  the  unloading  cycle.  Also,  in  the  finite  element  scheme,  the  crack 
surface  at  contact  was  fixed  against  further  displacement  in  the  plane  normal 
to  the  crack  surface,  whereas  the  cracks  were  allowed  to  interpenetrate  in  the 
BEM  modeling.  While  the  first  modeling  inconsistency  is  not  expected  to  have 
substantial  influence  on  the  results,  the  second  modeling  difference  may  be 
the  cause  of  the  discrepancies  in  the  results,  since  the  premature  closing  of 
crack  changes  the  crack  tip  locations.  The  crack  tip  singularity  which  drives 
the  reversed  plastic  strain  effectively  vanishes  at  the  initial  crack  tip,  as 
a  consequence,  the  development  of  reversed  plastic  strain  is  greatly 
restrained.  Nevertheless,  the  fact  that  different  crack  opening/closing 
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25.  Comparison  of  Crack  Opening/Closing  Stresses  for 
Constant  K  and  Constant  S  Cases — Plane  Strain 


behavior  for  different  loading  situation  (i.e.,  constant  S(max)  vs.  constant 
K(max))  using  the  same  approach  suggests  different  correlation  between  crack 
opening/closing  stress  and  maximum  stress  amplitude. 


Of  the  crack  cases  considered,  two  situations  corresponding  to 
crack  lengths  of  10  in.  and  0.1  in.  were  studied  in  detail.  The  long  crack 
case  corresponds  to  a  ratio  of  applied  maximum  stress  to  yield  stress  of  about 
0.71  and  the  short  case  corresponds  to  a  ratio  of  0.07.  Given  a  reasonable 
assumption  that  cyclic  crack  extension  is  associated  with  crack  tip 
deformation,  thereby,  with  the  crack  opening  displacements,  the  results 
reported  earlier  for  the  long  crack  plane  strain  constant  K  case  (Figure  13) 
predict  that  the  crack  responds  in  accordance  with  LEFM.  The  same  conclusion 
can  be  drawn  from  the  normalized  crack  opening  displacement  for  the  constant 
stress  case  plotted  in  Figure  26. 


The  short  crack  results  in  Figure  14  and  16  for  the  constant  K  case  and 
in  Figures  27  and  28  for  the  constant  stress  case  show  the  effect  of  plastic 
wake  on  COD  at  the  region  near  the  crack  tip  as  a  result  of  high,  unreversed 
crack  tip  stretch  plasticity.  The  plane  stress  results  for  both  long  (Figures 
15  for  constant  K,  and  Figure  29  for  constant  S)  and  short  (Figure  16  for 
constant  K,  and  Figure  28  for  constant  S)  cracks  show  enduring  plastic  wake 
for  both  cases.  The  higher  stretch  within  the  scale  of  the  cyclic  crack  tip 
zone  than  the  LEFM  prediction  may  provide  the  mechanics  for  increased  short 
crack  fatigue  crack  growth. 


On  the  other  hand,  the  measure  of  the  closure,  evaluated  in  terms  of 
crack  opening/closing  stress  levels,  does  not  seem  to  provide  the  explanation 
for  fatigue  crack  growth.  The  amount  of  plastic  wake  affects  crack  closure, 
but  the  closure  seems  not  to  be  the  cause  of  increased  fatigue  crack  damage. 
Yet,  the  amount  of  closure  is  an  indication  of  the  excess  crack  tip  damage  and 


P'one  Stress.  Long  Crack,  a='0.1 


0  25*S(max)  +  0.5*S(max) 
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thus  serve,  in  a  secondary  manner,  as  a  correlating  parameter  for  short  crack 
growth  rate  data.  It  is  hypothesized  that  stretch,  rather  than  closure  is  the 
concern  in  predicting  the  short  crack  growth  rates  observed  in  testing. 


3.0  THREE-DIMENSIONAL  FRACTURE  MECHANICS  MODELING 


3.1  Overview  of  LEFM  and  BEM  Analysis 

Material  deficiencies  in  the  form  of  pre-existing  flaws  initiate  cracks 
and  fractures  in  structures.  The  presence  of  a  crack  in  structural  element 
generally  induces  high  stress  concentration  at  the  crack  tip  and  thereby 
reduces  the  strength  of  the  structure.  Fracture  mechanics  provide 
satisfactory  means  for  the  characterization  of  this  local  crack  tip  stress 
fields  as  well  as  the  elastic  deformations  of  the  material  in  the  neighborhood 
of  the  crack.  In  LEFM,  the  inelastic  deformation  in  the  vicinity  of  the  crack 
tip  due  to  stress  concentrations  is  deemed  to  be  small  compared  to  the  size  of 
the  crack  and  other  characteristic  lengths. 

Elastic  modeling  of  crack  tip  behavior  makes  use  of  deformation  due  to 
three  primary  modes  of  loading  as  illustrated  in  Figure  30.  The  three  modes 
are:  the  opening  mode  (Mode  I)  due  to  normal  stress,  sliding  mode  (Mode  II) 
due  to  in-plane  shear  stress,  and  tearing  mode  (Mode  III),  due  to  out-of-plane 
shear. 

Consider  the  crack  problem  of  Figure  31,  representing  an  infinite  plate 
under  triaxial  stress.  The  stresses  and  displacements  for  traction  free 
cracks  may  be  given  as  an  infinite  series  in  r,  where  r  is  the  distance  from 
the  crack  tip.  For  the  antiplane  problem  (see  [6]),  the  near  crack  tip  field 
is  given  by: 
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f  3U  .  _ UJ— .  f  1  (5) 

°32  (2,r)1/2  cos(®/2r 


and 


2K 

II I f  r  >1/2 
u  ^2*' 


sin(e/2) 
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where  is  the  displacement,  »tj  is  the  stress  and  n  is  the  shear  modulus. 

The  Mode  III  stress  intensity  factor,  Kjh,  is  established  by  the  far  field 
boundary  conditions  and  is  a  function  of  applied  loading  and  the  geometry  of 
the  cracked  body.  The  Mode  III  stress  intensity  factor  is  defined  by: 

K  =  lim{(2*r)1/2a,.  }  (7) 

111  r*o  32  9=0 

For  the  plane  problem,  assuming  plane  strain  conditions, 


cos( 9/2) 


(2irr) 


1  -  sin(0/2)sin(39/2) 
sin( 9/2)cos(38/2) 

1  «■  sin(0/2)sin(39/2) 


Ul  Kr  r  1/2  \  cos(6/2)  [k  -  1  ♦  2sin  (0/2) ]  l 

u2  2u  2*  1  sin(8/2)  [k+1-2cos2(0/2) 1  ) 


where  Kj  is  the  Mode  I  stress  intensity  factor,  defined  by: 


Cr  *  lim  {(2*r)1/2o. 


k  =  3-4v 


The  corresponding  relationships  for  Mode  II  field  are: 
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where  Mode  II  stress  intensity  factor  Kji.  is  defined  by: 


r , -  ,1/2 

K  .  =  Urn  ((2*r) 

11  r*0 
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3.1.1  Boundary  Element  Method  for  Fracture  Mechanics  Analysis 

The  analytical  basis  of  the  method  is  the  transformation  of  the 
governing  equilibrium  equation  of  an  isotropic,  homogeneous,  elastic  element 
by  an  integral  identity,  using  Betti's  reciprocal  work  theorem.  The  identity 
for  the  displacement  at  a  point  P(x)  is  given  by  (e.g.,  [34]): 


CtJ(P)u  <P) 


l  T  <P,Q)u  (Q)ds(Q) 
S  J  J 


Ju. ,(P,Q)t  (Q)ds(Q)  (14) 
S  J 


where  t^(Q)  and  ui(Q)  are  the  boundary  values  of  traction  and  displacement, 
Tij(P.Q)  and  Ujj(P.Q)  are  tractions  and  displacements,  respectively,  in  xt 
directions  at  Q(x)  due  to  orthogonal  unit  loads  in  the  xj  directions  at  P(x) . 
The  discontinuity  term  Cjj  is  equal  to  1/2  for  smooth  boundary  points  and  can 
be  evaluated  indirectly  using  rigid  body  translation  as  described  by  Cruse 
[42],  for  non-smooth  boundary  points. 

The  utility  of  the  method  as  a  general  practical  solution  tool  is 
facilitated  by  two  approximations;  one  is  the  description  of  the  boundary  S  by 
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a  finite  number  of  surface  elements,  the  second  is  the  representation  of  the 
field  variables  (u^.t^)  and  geometry  by  known  interpolation  functions  within 
individual  elements.  In  the  present  analysis,  triangular  and  quadrilateral 
elements  are  used  for  surface  representation.  The  field  variables,  as  well  as 
geometry,  are  represented  by  isoparametric  quadratic  interpolation  functions. 
Numerical  evaluation  of  discretized  integrals,  as  well  as  the  use  of  special 
crack  tip  elements  employed  in  the  current  analysis,  are  described  in  a  later 
section. 

3.1.2  BEM  Modeling  of  Cracked  Bodies 


The  numerical  solution  of  Eq.  (14)  is  straight-forward  for  a 
general  three-dimensional  stress  analysis.  However,  the  presence  of  two 
coplanar  surfaces  preclude  the  use  of  the  method  for  general  solution  of 
cracked  bodies.  Therefore,  several  different  modeling  strategies  have  been 
employed  for  three-dimensional  cracked  bodies,  as  illustrated  in  Figure  32. 

The  first  approach  is  to  model  the  crack  as  an  open  notch  as  reported  by  Cruse 
(41).  The  major  deficiency  of  this  modeling  approach  is  that  the  results 
become  inaccurate  when  the  surfaces  are  modeled  too  far  apart,  however,  the 
system  equation  becomes  badly  conditioned  when  the  surfaces  are  modeled  too 
close  together.  One  form  of  avoiding  this  difficulty  is  the  dislocation  or 
traction  BIE  modeling  approach  as  developed  in  different  forms  by  Cruse  [42], 
Guidera  and  Lardner  [45],  Bui  [38]  and  Weaver  [49],  The  singular  nature  of 
the  integrals  is  this  method  poses  difficulty  in  the  numerical  implementation 
as  reported  by  these  authors.  Significant  improvements  have  been  reported 
recently  by  Polch,  et  al.  [48]  and  others.  However,  further  research  is 
required  before  the  full  potential  of  this  utility  can  be  realized. 

Modeling  of  symmetric  cracks  is  rather  straight-forward  since  only 
a  symmetric  part  of  the  body  that  contains  one  crack  surface  needs  to  be 
modeled.  Earlier  application  of  this  modeling  has  been  reported  by  Cruse  and 
VanBuren  [40],  Cruse  [42],  and  Cruse  and  Meyers  [44].  These  results  have  been 
subsequently  improved  by  using  isoparametric  interpolation  functions  by  Cruse 
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Figure  32.  BEM  Modeling  Strategies  for  Three-Dimensional  Cracks 


and  Wilson  [43] .  The  accuracy  of  the  results  has  been  further  improved  by  the 
use  of  special  crack  tip  elements  that  accurately  model  the  crack  tip  field. 
Improved  results  using  higher  order  interpolation  functions  for  geometry  and 
field  variables  and  special  crack  tip  elements  are  also  reported  by  Cruse  and 
Wilson  [43].  A  detailed  description  of  special  crack  tip  elements  is  given  in 
the  following  section. 

A  modeling  strategy  applicable  for  a  general  three-dimensional 
non-symmetric  crack  is  the  sub-region  model.  In  this  approach,  the  body  is 
substructured  into  separate  regions  through  the  crack  plane  such  that  each 
crack  surface  is  in  a  different  region.  The  overall  solution  is  obtained  by 
satisfying  compatibility  and  continuity  conditions  at  the  interface  of  the 
regions  except  along  the  crack  surface.  Numerical  results  using  this  strategy 
for  two-dimensional  structures  have  been  reported  by  Blandford,  et  al.  [37]. 
Results  for  three-dimensional  bodies  are  reported  in  Section  3.2. 

3.1.3  Use  of  Singular  Elements 

The  accuracy  of  the  numerical  computations  is  enhanced  by  proper 
representation  of  the  field  variables  in  the  vicinity  of  the  crack  tip.  It  is 
well-known  that  the  crack  tip  opening  displacement  varies  with  square  of  the 
distance  (r)  from  the  crack  tip,  whereas,  the  stresses  produce  l//r 
singularity.  The  crack  tip  elements  are  modified  such  that  they  capture  these 
variaM  nns . 

It  has  been  observed  by  Barsoum  [36]  that  the  placement  of 
midpoint  nodes  of  the  element  sides  emanating  from  the  crack  tip  at  quarter 
points  leads  to  the  required  displacement  and  stress  variation  in  the  finite 
element  method.  By  using  quarter  point  (QP)  elements  (Figure  33)  on  both 
sides  of  the  crack  in  BEM,  the  displacements  and  tractions  are  made  to  vary  in 
physical  space  as: 

u(r)  ,  ,  .  /icx 

{.  }  =  A  ♦A2/(r/l)*A3(r/l)  0-5) 

t(r) 


JfX 


where  r  is  the  normal  distance  from  the  crack  tip,  1  is  the  element  length  in 
that  direction,  and  A^,  A2,  A3  are  functions  of  nodal  displacements  or 
tractions.  However,  unlike  the  finite  element  method  in  which  the  tractions 
(or  stresses)  are  derived  from  the  spatial  derivative  of  displacements,  the 
displacements  and  tractions  are  independently  approximated  in  BEM.  The  use  of 
quarter  point  elements  on  both  sides  of  the  crack  front,  thus,  will  give  the 
same  variation  for  both  displacements  and  tractions  and,  therefore,  Eq.  (15) 
does  not  give  the  required  singularity  for  crack  tip  tractions.  The  simplest 
way  to  obtain  this  singularity  appears  through  the  use  of  singular  shape 
functions.  However,  the  numerical  integration  scheme  employed  in  the  current 
analysis  requires  subsegmentation  of  elements,  as  described  in  the  following 
section,  which  makes  the  implementation  difficult.  Instead,  the  traction 
singularity  is  contrived  by  the  multiplication  of  shape  function  by  a 
non-dimensional  parameter  /(t/r)  as: 

t(r)  =  /( l/r ) -t  (16) 


where  t  is  the  normal  traction  defined  by  equation  (15).  The  variation  of 
traction  in  these  traction  singular  (TS)  elements  is  then  given  by: 


E(r)  !  7777TT  *  a2- 


The  use  of  these  special  elements  improve  the  accuracy  substantially  as  seen 
by  the  numerical  results  reported  herein. 
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3.1.4  Stress  Intensity  Factor  Evaluations 

The  stress  intensity  factors  may  be  computed  using  displacements 
or  traction  BEM  solution  from  Eqs.  5-14.  Mode  I  and  Mode  II  stress  intensity 
factors  are  computed  from: 


"I  or  II  = 


u  ,  / 2t . 

4(  1-v)  r  ;uI  or 


II 


(r) 


(18) 


where  uj  and  ujj  are  decoupled  displacements  along  Mode  I  and  Mode  II 
directions  at  the  quarter  point,  and  r  is  the  distance  of  the  quarter  point 
from  the  crack  front.  Mode  III  stress  intensity  factors  may  be  computed 
similarly  from: 


KIII  *  i  '  <r,u 1„ <r>  (l9) 

Alternatively,  the  stress  intensity  factors  can  be  computed  from 
decoupled  tractions  at  the  crack  front.  As  an  example,  the  Mode  III  stress 
intensity  factor  may  be  evaluated  from: 

Km  =  /(2irl)tm  (20) 

where  f,„  is  the  decoupled  nominal  traction  defined  by  Eq.  (16).  However,  in 
the  present  analysis,  the  stress  intensity  factors  are  computed  using  quarter 
point  displacements  (Eqs.  (18)  and  (19)),  since  the  values  using  crack  front 
tractions  (Eq.  (20))  generally  over-estimated  the  stress  intensity  factors  by 
about  10%.  It  is  believed  that  one  of  the  reasons  for  the  discrepancy  may  be 
that  higher  accuracy  is  needed  for  the  integration  of  traction  singular 
element  than  the  one  used  in  the  present  analysis.  Another  reason  may  be 
that,  though  the  first  term  in  Eq.  (17)  provides  the  required  singularity,  the 
higher  order  terms  of  standard  quadratic  element  expression  (15),  through 
traction  singular  modification,  are  the  suitable  ones  for  the  representation 


of  traction  variation.  However,  the  use  of  traction  singular  modification 
with  quarter  point  elements  improved  the  accuracy  of  the  displacement  based 
stress  intensity  factors  and  are,  therefore,  used  in  the  current  analysis. 

3.1.5  Evaluation  of  Discretized  Boundary  Integrals 

In  the  present  analysis,  the  geometry  as  well  as  field  variables 
are  represented  by  isoparametric  quadratic  shape  functions.  The  use  of  higher 
order  interpolation  functions,  in  general,  preclude  the  use  of  analytical 
integration  and,  therefore,  numerical  quadrature  is  used  in  the  current 
development.  Non-singular  Kernel  function -shape  function  products,  in 
principle,  can  be  directly  approximated  by  the  application  of  Gauss -Legendre 
quadrature  formula.  However,  to  maintain  a  certain  level  of  accuracy,  element 
subdivision  may  be  required.  Following  Lachat  and  Watson  [46] ,  the  minimum 
element  side  length  for  a  given  error  tolerance  and  quadrature  order  is 
determined  from  error  analysis.  The  element  is  then  subdivided  to  satisfy 
this  requirement. 

The  non- singular  integration  is  performed  through  a  polar 
coordinate  transformation  which  eliminates  the  singularity.  To  accomplish 
this,  the  element  is  subdivided  through  the  singular  point  and  the  polar 
coordinate  system  is  constructed  through  the  singular  apex  as  described  by 
Banerjee  and  Raveendra  [35]. 

To  achieve  the  required  crack  tip  field  variation,  considerable 
mesh  refinement  at  crack  tip  is  necessary.  Because  of  this  non-gradual  mesh 
refinement,  the  subdivision  scheme  employed  for  non- singular  integration  does 
not  always  work  efficiently.  In  order  to  avoid  a  large  number  of 
subdivisions,  a  modified  singular  integration  scheme  is  used  for  these 
non- s ingular  integrations . 
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3 . 2  Numerical  Modeling  of  Elastic  Crack 

The  solution  procedure  outlined  in  the  previous  section  has  been  applied 
for  the  solution  of  many  symmetric  and  non-symmetric  cracked  specimens.  The 
numerical  solutions  obtained  by  using  the  current  computer  code  CRX3D  are 
validated  against  known  analytical  solutions,  in  some  cases,  and  are  compared 
to  finite  element  solutions  in  other  cases. 

3.2.1  Circular  Crack 

Figure  34  shows  the  BEM  map  used  for  the  analysis  of  a  buried 
circular  crack  problem.  Due  to  symmetry,  only  one-eighth  of  the  body  was 
modeled.  In  addition,  the  symmetric  planes  of  the  body  were  not  modeled.  The 
crack  opening  displacements  using  standard  quadratic  element  and  modified 
crack  tip  elements  are  compared  to  analytical  results  in  Figure  35  which 
indicate  that  the  results  are  improved  by  the  use  of  modified  crack  tip 
elements.  The  analytical  solution  plotted  in  the  figure  is  for  an  infinite 
body.  To  further  assess  the  accuracy,  the  Mode  I  stress  intensity  factors 
were  evaluated  for  standard  and  modified  crack- tip  element  models.  A 
comparison  with  the  empirical  value  that  takes  the  finite  dimension  of  body 
into  account  demonstrated  that  the  stress  intensity  factors  using  quarter 
point  element  (QP)  is  approximately  5%  in  error,  compared  to  standard  element 
results  which  is  approximately  10%  in  error.  Further  enhancement  was  obtained 
by  using  traction  singular  (TS)  modification  to  the  quarter  point  element, 
which  improved  the  stress  intensity  factor  value  to  within  1%  of  the  predicted 
value . 

3.2.2  Elliptical  Surface  Crack 

Figure  36  shows  the  BEM  map  for  one -fourth  of  a  finite  cracked 
plate.  Whereas  the  symmetric  Y-Z  plane  was  not  modeled,  the  X-Z  plane  was 
modeled  to  allow  the  solution  of  elliptic  buried  and  semi-elliptic  surface 
cracks  to  be  modeled  by  a  change  in  boundary  condition.  Figure  37  shows  the 
normalized  Mode  I  stress  intensity  factor  for  buried  crack  using  modified 
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Figure  34.  Cricular  Crack — BEM  Map 
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Figure  37.  Mode  I  Stress  Intensity  Factor  for  Buried  Elliptical  Crack 
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crack  tip  elements  compared  to  the  finite  body  solution.  The  figure  indicates 
that  the  accuracy  of  the  solution  was  again  improved  by  the  use  of  the 
traction  singular  model. 

The  semi-elliptic  surface  crack  model  was  subjected  to  both 
tensile  and  bending  loads.  Figures  38  and  39  show  good  agreement  between 
CRX3D  results  and  the  empirical  Mode  I  stress  intensity  factors  provided  by 
Newman  and  Raju  [47]  for  both  loading  cases. 

3.2.3  Inclined  Circular  Crack 

Due  to  the  symmetry,  only  one  surface  of  the  crack  was  modeled  in 
the  previous  examples.  However,  as  explained  in  the  previous  section,  the 
non-symmetric  cracks  were  modeled  using  subregion  modeling  strategy.  CRX3D 
code  was  validated  for  non-symmetric  crack  cases  by  analyzing  a  circular  crack 
which  is  inclined  45%  to  the  loading  direction.  Figure  40  shows  the  BEM  model 
of  one -half  of  the  cracked  body.  Stress  intensity  factors  for  all  three  modes 
compared  well  to  the  infinite  body  analytical  solution,  as  shown  in  Figure  41. 


3.2.4 


The  final  example  solved  is  a  T-Joint  section  (Figure  42(a))  that 
comprises  a  part  elliptic  surface  flaw.  A  BEM  model  of  one-half  of  the  body 
is  shown  in  Figure  42(b)).  Again,  a  substructure  modeling  strategy  was  used. 
The  body  was  modeled  into  four  subregions  to  further  improve  the  accuracy. 
Figure  44  shows  the  Mode  I  stress  intensity  factor  normalized  with  respect  to 
a  two-dimensional  plane  strain  model  result  (Figure  43).  The  behavior 
indicates  that  while  the  stress  intensity  factor  is  higher  at  the  free-surface 
of  the  body,  the  solution  at  the  mid- surface  is  comparable  to  the 
two-dimensional  solution  as  expected.  Mode  II  and  Mode  III  stress  intensity 
factors  are  also  computed  for  the  model  and  are  shown  in  Figure  45. 
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Figure  38.  Mode  I  Stress  Intensity  Factor  for 
Seai-Elliptic  Surface  Crack-Tensile  Load 
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Figure  39.  Mode  I  Stress  Intensity  Factor  for 
Seai-Elliptic  Surface  Crack- -Bending  Load 
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Figure  45.  Mode  II  and  Mode  III  Stress  Intensity  Factors  for  T-Joint  Model 


As  seen  in  the  previous  section,  unlike  the  two-dimensional  approach  in 
which  the  use  of  special  fundamental  solution  that  satisfies  the  stress  free 
condition  at  the  crack  surface  eliminates  the  crack  surface  modeling,  crack 
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surfaces  are  modeled  in  the  three-dimensional  approach.  In  modeling  the 
three-dimensional  crack  problem,  the  surface  of  the  body  is  discretized  into 
quadrilateral  and  triangular  elements  over  which  the  variables  are 
approximated  by  quadratic  interpolation  functions.  The  interior  plastic 
strains  are  modeled  using  rectilinear  cells  over  which  the  plastic  strains  are 
assumed  constant.  Due  to  the  use  of  higher  order  interpolation  functions,  the 
surface  integrals  in  BIE,  which  cannot  be  evaluated  analytically,  are 
evaluated  numerically.  The  volume  integrals,  however,  are  evaluated 
analytically.  The  elastoplastic  boundary  element  solution  approach  involves 
iterative  solution  of  the  displacement  as  well  as  the  displacement  gradient 
forms  of  the  boundary  integral  equation,  as  described  in  detail  in  the  first 
annual  report.  The  boundary  integral  equation  (BIE)  corresponds  to  the 
displacement  gradient  is  deduced  from  the  derivative  of  the  displacement 
equation,  thus,  the  singularity  of  the  integrands  are  increased.  Generally, 
the  evaluation  of  these  hyper -singular  integrals  are  numerically  challenging 
for  both  singular  and  near  singular  cases.  Especially  for  crack  problems,  one 
must  have  small  domain  elements  for  proper  plastic  strain  modeling.  The  need 
to  compute  the  near  singular  displacement  gradient  integrals  at  the  centroid 
of  the  domain  cells  which  are  close  to  the  crack  surface  imposes  a  heavy 
burden  on  the  analyst.  While  adaptive  transformation  techniques  are  used  to 
ease  the  numerical  burden  and  improve  the  accuracy  of  these  integrals  to  some 
extent,  the  magnitude  of  the  problems  is  still  substantial.  Due  to  the 
constraint  imposed  by  the  need  to  accurately  evaluate  the  hyper -singular 
integrals  and  to  limit  the  overall  computing  size  a  somewhat  coarse  mesh  was 
used  in  the  current  analysis. 
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The  model  investigated  in  the  current  study  is  that  of  a  center,  through 
crack  of  length  0.2  in.  in  a  square  plate  of  length  20  in,  and  thickness  2  in. 
Due  to  symmetry,  only  one  quadrant  of  the  cracked  body  is  modeled.  This  model 
corresponds  to  the  two-dimensional  short  crack  case  investigated  in  Section 
2.0. 

Adequate  modeling  of  crack  tip  requires  small  elements  near  the  crack  tip 

region  of  a  crack  that  is  small  compared  to  the  overall  geometry  of  the  body. 

Since  a  single  region  boundary  element  modeling  approach  will  affect  the 
numerical  stability  of  the  solution  adversely,  the  crack  was  modeled  into  two 
regions ,  an  inner  region  around  the  crack  as  shown  in  Figure  46 ,  and  an  outer 
region  as  shown  in  Figure  47.  The  interior  plastic  strains  were  modeled  using 

the  cells  shown  in  Figure  48.  The  material  properties  used  are  the  same 

values  used  for  the  two-dimensional  case.  Only  the  constant  amplitude  stress 
loading  is  considered.  Loading  consisted  of  incremental  loading  up  to  a 
maximum  value  corresponding  to  20  ksi./in  and  subsequent  unloading  to  zero 
value.  Crack  is  simulated  by  imposing  zero  boundary  condition  normal  to  the 
crack  plane  along  the  plane  that  contains  the  crack  except  at  the  crack 
surface.  Crack  extension  is  simulated  by  changing  the  boundary  condition  to 
the  appropriate  new  crack  locations. 

3.4  Numerical  Results  for  Elastoplastic  Cracks 

During  the  first  loading  cycle,  the  crack  opening  displacements  were 
monitored  along  the  crack  surface.  The  crack  opening  displacement  normalized 
with  respect  to  the  elastic  value  at  the  center  of  the  crack  is  plotted  in 
Figure  49.  The  crack  opening  displacement  at  K  -  5  ksi./in  during  the  loading 
cycle  is  essentially  elastic  due  to  the  finite  size  of  the  smallest  domain 
element  used  for  the  plastic  strain  modeling.  At  the  maximum  load  of  K  -  20 
ksi./in,  the  plasticity  causes  the  displacement  to  be  higher  than  the  elastic 
value  and  this  is  seen  for  all  crack  points.  The  ratio  of  the  crack  opening 
displacement  computed  during  the  unloading  cycle  is  higher  since  only  a 
portion  of  the  plastic  stretch  is  reversed  during  unloading.  The  same  plot  at 
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Figure  49.  Normalized  COD  at  Initial  Crack  Location — y/a=10 
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the  free  surface  shown  in  Figure  50  indicates  similar  behavior.  However,  as 
seen  in  Figure  51,  the  crack  opening  displacements  increase  progressively  from 
the  center  of  the  crack  (y/a-10)  to  the  free  surface  (y/a-0) . 

The  crack  was  extended  by  finite  amounts  of  0.125a  and  the  loading  cycle 
was  repeated  at  each  of  the  new  crack  length.  Equilibrium  at  the  new  crack 
length  was  achieved  by  solving  an  elastic  problem  with  previous  plastic 
strain,  by  a  single  iterative  cycle.  Figure  52  plots  the  crack  opening 
displacement  normalized  to  the  elastic  values  at  the  center  crack  during  the 
unloading  cycle  after  extending  the  crack  by  0.25a.  The  situation  at  the 
center  of  the  crack  approaches  plane  strain  conditions,  however,  unlike  the 
plane  strain  results  reported  earlier,  these  results  indicate  no  closure.  The 
same  plot  at  the  free  surface,  Figure  53,  where  the  condition  corresponds  to 
plane  stress  situation,  also  indicates  no  closure.  These  results  are 
inconsistent  with  the  expectations.  The  cells  yielded  at  the  maximum  and 
minimum  loads  are  plotted  in  Figure  54.  The  stresses  at  no  load  during 
unloading  is  compressive,  therefore,  it  is  reasonable  to  expect  closure. 

The  crack  opening  displacement  nearest  to  the  crack  tip  during  the 
initial  and  final  loading  cycles  are  plotted  in  Figure  55  at  the  center  of  the 
crack.  The  displacements  are  normalized  to  the  elastic  value  at  the  maximum 
load.  The  results  show  that  when  extended  at  zero  load  the  crack  barely 
closes  even  though  compressive  stresses  are  present.  During  the  subsequent 
loading,  the  combination  of  reduced  stiffness  and  small  size  of  reversed 
plasticity  compared  to  the  large  amount  of  forward  plasticity  seems  to  cause 
the  crack  to  remain  open  at  no  load.  Similar  results  are  obtained  for  the 
crack  tip  node  at  the  free  surface  (see  Figure  56).  Since  the  size  of  the 
plastic  cells  used  in  the  three  dimensional  model  case  is  much  larger  than 
the  two-dimensional  model,  a  comparable  two-dimensional  model  was  created  to 
analyze  the  influence  of  mesh.  However,  the  two-dimensional  results  for  the 
plane  stress  and  plane  strain  cases,  shown  respectively  in  Figure  57  and  58, 
indicate  closure. 
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Figure  50.  Normalized  COD  at  Initial  Crack  Location — y/a=0 
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Figure  51.  Variation  of  COD  Through  Thickness 
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Pigure  52.  Normalized  COD  After  Crack  Extension —  y/a*  10 


In  modeling  the  two-dimensional  crack,  the  elastic  solution  to  the  crack 
problem  is  'exact',  and  the  only  approximation  is  the  domain  discretization 
used  for  the  plastic  strain.  In  the  three-dimensional  case,  however,  both  the 
elastic  crack  problem  and  the  plastic  strains  are  approximated.  Further,  the 
discretized  integrals  in  the  two-dimensional  case  were  evaluated  analytically, 
whereas,  the  three-dimensional  integrals,  except  the  volume  integrals,  were 
evaluated  using  approximate  numerical  procedures  due  to  the  use  of  higher 
order  interpolation  functions.  While  the  use  of  higher  interpolation 
functions  for  the  surface  modeling  is  known  to  give  more  accurate  results,  as 
evidenced  in  Section  3.2,  the  need  to  compute  the  burdensome  near  singular 
hyper- singular  integrals  accurately  for  displacement  gradient  appears  to  be 
the  cause  for  the  inconsistencies  of  the  results.  It  is  believed  that  the 
discrepancies  between  the  two-  and  three-dimensional  results  are  mainly  a 
numerical  problem,  rather  than  a  difference  in  the  mechanics  of  the  problems. 
Further  research  is  required  to  refine  the  three-dimensional  modeling  problem 
to  complete  the  analysis. 


4.0  PROGRAM  ACCOMPLISHMENTS  AND  CONCLUSIONS 


The  principal  accomplishment  of  the  project  is  the  development  of  a 
boundary  element  method  analysis  to  investigate  the  effect  of  plasticity  on 
fatigue  crack  propagation  of  two-  and  three-dimensional  cracks.  The 
two-dimensional  approach,  based  on  the  use  of  a  special  fundamental  solution 
which  explicitly  accounts  for  the  presence  of  stress -free  crack,  was  developed 
in  a  previously  sponsored  AFOSR  contract.  The  approach  is  extended  for  the 
investigation  of  crack  tip  behavior  during  crack  propagation  under  cyclic 
loading.  For  the  three-dimensional  case,  a  previously  developed  computer 
program  was  modified  to  model  fracture  mechanics  problems,  this  includes  the 
use  of  substructuring  for  non- symmetric  cracks  and  the  use  of  special 
crack- tip  elements.  The  elastoplastic  analysis  uses  volume  discretization  and 
the  volume  contribution  of  the  plastic  strain  is  computed  using  an  analytical 
solution  of  a  rectilinear  cell  that  contains  plasticity  embedded  in  an 
infinite  body.  The  major  findings  of  the  project  are  as  follows: 

(1)  The  plastic  zone/residual  plastic  zone  size  for  stationary  and 
propagating  cracks  is  estimated  by  a  new  indirect  method.  This  is,  by 
computing  the  elastic  stress  intensity  factor,  due  to  crack  extension 
into  prior  plasticity,  the  influence  region  of  residual  plasticity  is 
estimated.  This  gives  an  unambiguous  estimate  of  the  plastic  zone  and 
compares  favorably  with  the  Dugdale  model  prediction  for  small  plastic 
zones . 

(2)  It  is  established  that  the  effect  of  the  plastic  wake  on  the  stress 
intensity  factor  for  crack  opening  (closure)  is  identical  to  the  effect 
of  the  residual  stress  of  crack  on  the  retardation  of  the  stress 
intensity  factor.  This  is,  crack  retardation  and  closure  are  seen  as 
identical  manifestation  of  the  same  plasticity  process  for  both 
stationary  and  propagating  cracks.  Further,  the  BIE  formulation  gives 
insight  to  the  mathematical  equivalence  of  these  two  effects. 


(3)  It  is  proposed  plastic  stretch  rather  than  crack  closure  is  the  cause  of 
increased  fatigue  crack  damage  of  small  cracks.  The  loss  of  constraint 
going  from  long  crack  to  short  crack  which  is  further  diminished  going 
from  plane  strain  to  plane  stress  causes  increase  in  plastic  stretch  and, 
therefore,  invalidity  of  LEFM.  Nevertheless,  the  amount  of  closure  is  an 
indication  of  the  excess  crack  tip  plastic  damage  and,  thus,  may  serve  in 
a  secondary  manner  as  a  correlating  parameter  for  short  crack  fatigue 
growth  rate  data. 

(4)  Three-dimensional  elastoplastic  problem  imposes  high  requirement  on 
finding  efficient  means  for  the  accurate  evaluation  of  integrals  and  also 
the  size  of  the  problem  makes  it  computationally  burdensome  for  any 
numerical  solution  procedure. 

(5)  The  results  establish  BEM  as  a  very  attractive  solution  tool  for  the 
elastic  and  elastoplastic  fracture  mechanics  analysis,  especially  for  the 
two-dimensional  case.  Since  only  the  plastic  strain  around  the  crack 
tip,  rather  than  the  entire  field,  as  is  the  case  in  finite  element 
method,  is  approximated  in  the  boundary  element  method,  the  results  are, 
generally,  more  accurate.  Even  though  the  domain  modeling  required  for 
the  elastoplastic  solution  may  appear  to  make  the  method  inefficient,  the 
fact  is  that  only  the  yielded  region  around  the  crack  needs  domain 
modeling  and,  therefore,  compared  to  the  overall  problem  the  domain 
modeling  is  not  substantial. 

(6)  Finally,  the  research  results  have  been  presented  in  a  variety  of 
technical  symposia  and  will  result  in  archival  articles.  The  list  of 
articles  and  presentations  resulted  form  the  research  is  as  follows: 

1.  "A  Comparison  of  Long  and  Short  Crack  Elastoplastic  Response  Using  the 


Boundary  Element  Method,"  T.A.  Cruse  and  S.T.  Raveendra, 
Fracture  Mechanics,  to  be  published. 
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2.  "BEM  Analysis  of  Problems  of  Fracture  Mechanics,"  S.T.  Raveendra  and  T.A. 
Cruse,  (Editors  P.K.  Banerjee  and  R.B.  Wilson)  Developments  in 
Boundary  Element  Methods  -5.  to  be  published. 

3.  "Application  of  the  Boundary  Element  Method  to  Inelastic  Fracture 
Mechanics,"  T.A.  Cruse  and  S.T.  Raveendra,  (Proceedings  of  the 
International  Conference  on  Computational  Engineering  Science,  Atlanta, 
Georgia,  April  10-14,  1988). 

4.  "Some  Elastoplastic  Fracture  Mechanics  Results  Using  Boundary  Integral 
Equations,"  presented  at  First  World  Congress  on  Computational  Mechanics, 
University  of  Texas  at  Austin,  September  22-26,  1986. 

5.  "Boundary  Element  Fracture  Mechanics  Modeling,"  presented  at  the  Fourth 
International  Conference  on  Numerial  Methods  in  Fracture  Mechanics,  San 
Antonio,  Texas,  March  23-27,  1987. 

6.  "Application  of  the  Boundary  Element  Method  to  Nonlinear  Stress 
Analysis,"  presented  at  ASCE  Engineering  Mechanics  Division  Specialty 
Conference,  State  University  of  New  York  at  Buffalo,  New  York,  May  20-22, 
1987. 

7.  "Boundary  Element  Fracture  Mechanics  Modelling,"  presented  at  the  IUTAM 
Symposium  on  Advanced  Boundary  Element  Methods:  Applications  in  Solid 
and  Fluid  Mechanics,  San  Antonio,  TX,  April  13-16,  1987. 
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